"""
2.	梅尔频率倒谱系数（MFCC：Mel Frequency Cepstral Coefficents）是一种在自动语音和说话人识别中广泛使用的特征，
它是由Davis和Mermelstein在1980年提出，从那时起，在语音识别领域，MFCC得到广泛应用。MFCC的提取过程如下图所示，
按照题目要求完成相应操作（44分）
"""
# ①	导入必要的语音处理工具包，导入并显示原始语音波形
import scipy.io.wavfile
import scipy.fftpack
import numpy as np
import os
import sys
import matplotlib.pyplot as plt
import seaborn as sns

BASE_DIR, FILENAME = os.path.split(__file__)
path = 'data/zsn-stop1.wav'
AUDIO_PATH = os.path.join(BASE_DIR, path)

sr, signal = scipy.io.wavfile.read(AUDIO_PATH)
if len(signal.shape) >= 2:
    signal = signal[:, 0]

spr = 2
spc = 1
spn = 0
plt.figure(figsize=[6, 6])
spn += 1
plt.subplot(spr, spc, spn)
plt.title('audio')
plt.plot(signal)
plt.grid()

# ②	图示预加重语音波形（预加重系数：0.97）
signal = np.append(signal[0:1], signal[1:] - 0.97 * signal[:-1])

spn += 1
plt.subplot(spr, spc, spn)
plt.title('audio emphasized')
plt.plot(signal)
plt.grid()
plt.show()

# ③	进行时域分帧加窗操作，并打印输出数据维度
signal_len = len(signal)
frame_len = int(np.ceil(sr / 40))
frame_stride = int(np.ceil(sr / 100))
n_frames = int(np.ceil((np.absolute(signal_len - frame_len) + 1) / frame_stride))
padded_len = n_frames * frame_stride + frame_len
print('signal_len', signal_len)
print('frame_len', frame_len)
print('frame_stride', frame_stride)
print('n_frames', n_frames)
print('padded_len', padded_len)
signal = np.append(signal, np.zeros(padded_len - signal_len, dtype=signal.dtype))

idx_col = np.arange(frame_len).reshape(1, -1)
idx_row = np.arange(0, frame_stride * n_frames, frame_stride).reshape(-1, 1)
idx_frame = idx_col + idx_row
frames = signal[idx_frame]

win = np.hamming(frame_len)
frames *= win
print('frames', frames.shape)

# ④	将短时语音进行快速傅里叶变换（FFT），转为频域分析，打印输出数据维度
NFFT = 512
rfft = np.fft.rfft(frames, n=NFFT)
rfft = np.absolute(rfft)
print('rfft', rfft.shape)
_, n_spec = rfft.shape
print('n_spec', n_spec)

spr = 1
spc = 2
spn = 0
plt.figure(figsize=[12, 6])
spn += 1
plt.subplot(spr, spc, spn)
plt.title('rfft')
sns.heatmap(rfft)

# ⑤	将数据幅度谱转化为功率谱，并显示功率谱2维图像
power_spec = rfft ** 2 / NFFT
print('power spectrum', power_spec.shape)
spn += 1
plt.subplot(spr, spc, spn)
plt.title('Power spectrum')
sns.heatmap(power_spec)
plt.show()

# ⑥	设置40个Mel三角滤波器组，进行Mel频率等间隔划分
n_filters = 40


def hz2mel(hz):
    return 2595 * np.log10(1 + hz / 700)


def mel2hz(mel):
    return (10 ** (mel / 2595) - 1) * 700


low_hz = 0
high_hz = sr / 2
low_mel = 0
high_mel = hz2mel(high_hz)
mel_arr = np.linspace(low_mel, high_mel, n_filters + 2)
hz_arr = mel2hz(mel_arr)
bins = np.floor(hz_arr / high_hz * n_spec)

spr = 2
spc = 1
spn = 0
plt.figure(figsize=[8, 8])
spn += 1
plt.subplot(spr, spc, spn)
plt.title('Hz')
y1y2 = np.tile([[100], [-100]], len(bins))
x1x1 = [bins, bins]
plt.plot(x1x1, y1y2)

# ⑦	根据Mel频率与普通频率之间的关系，确定40个Mel三角滤波器的频率范围，并显示Mel滤波器组2维图像
filters = np.zeros((n_filters, n_spec), dtype=np.float64)
print('filters', filters.shape)
for i in range(1, n_filters + 1):
    m_minus = int(bins[i - 1])
    m = int(bins[i])
    m_plus = int(bins[i + 1])
    for k in range(m_minus, m):
        filters[i - 1, k] = (k - bins[i - 1]) / (bins[i] - bins[i - 1])
    for k in range(m, m_plus):
        filters[i - 1, k] = (bins[i + 1] - k) / (bins[i + 1] - bins[i])
spn += 1
plt.subplot(spr, spc, spn)
plt.title('Filters')
for f in filters:
    plt.plot(f)
plt.show()

# ⑧	功率谱经过Mel三角滤波器组运算后，打印输出数据维度，并显示2维图像
filtered_power = np.dot(power_spec, filters.T)
print('filtered_power', filtered_power.shape)
spr = 1
spc = 3
spn = 0
plt.figure(figsize=[18, 6])
spn += 1
plt.subplot(spr, spc, spn)
plt.title('Filtered power')
sns.heatmap(filtered_power)

# ⑨	对数据谱进行对数运算获得数据倒谱
filtered_power = np.where(filtered_power <= 0., np.finfo(float).eps, filtered_power)
filtered_power = np.log10(filtered_power)
spn += 1
plt.subplot(spr, spc, spn)
plt.title('Logged')
sns.heatmap(filtered_power)

# ⑩	数据倒谱经过离散余弦变换（DCT）
dct = scipy.fftpack.dct(filtered_power, type=2, norm='ortho')
dct = dct[:, 1:13+1]

# 11	获得数据MFCC特征参数，并显示2维图像
spn += 1
plt.subplot(spr, spc, spn)
plt.title('MFCC')
sns.heatmap(dct)
plt.show()
